

# Setup ------------------------------------------------------------------------
  # Options
  Sys.setenv(TZ = "US/Pacific")
  options(stringsAsFactors = F)
  options(bitmapType = "cairo")
  # Packages
  library(pacman)
  p_load(
    lfe, prism, raster, data.table, lubridate, stringr, readr,
    ggplot2, ggforce, showtext, Cairo, rmapshaper, ggthemes, viridis,
    magrittr
  )
  try(ymd(today()))
  # Directories
  dir_project <- "/Users/edwardarubin/Dropbox/Research/MyProjects/NaturalGas/"
  dir_csv     <- paste0(dir_project, "DataCsv/")
  dir_rds     <- paste0(dir_project, "DataR/")
  dir_spatial <- paste0(dir_project, "DataSpatial/")
  dir_fig     <- paste0(dir_project, "Text/SubmissionJAERE/Figures/")
  dir_figures <- paste0(dir_project, "Text/SubmissionJAERE/Figures/")
  # Define colors
  very_light_grey = "#E0E0E0"
  light_grey      = "#BDBDBD"
  mid_grey        = "#9E9E9E"
  dark_grey       = "#616161"
  red_pink        = "#E40045"
  # purple          = "#310769"
  purple          = "#6A1B9A"
  light_purple    = "#9E17A0"
  orange          = "#FFA000"
  # Load fonts
  showtext_auto()
  font_add(family = 'Charter', regular = '/System/Library/Fonts/Supplemental/Charter.ttc')
  # My ggplot themes
  theme_beamer <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "Roboto", color = "black",
      size = 16),
    title = element_text(size = 18),
    legend.text = element_text(size = 18),
    legend.key.size = unit(2.5, "cm"))
  theme_paper <- theme_bw() + theme(
    panel.border = element_blank(),
    axis.ticks = element_blank(),
    legend.position = 'bottom',
    legend.key = element_blank(),
    panel.grid = element_blank(),
    text = element_text(family = "Charter", color = "black",
      size = 11),
    axis.text = element_text(color = "black"),
    title = element_text(size = 12),
    legend.text = element_text(size = 11),
    legend.key.size = unit(1.5, "cm"))
  # My own save functions
  save_beamer <- function(gg_tmp, name, width = 16 * 0.7, height = 10 * 0.7,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_beamer + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      # extrafont::embed_fonts(paste0(dir_figures, name, ".pdf"))
  }
  save_paper <- function(gg_tmp, name, width = 6.4, height = 4,
    themes = NULL) {
      ggsave(filename = paste0(name, ".pdf"),
        path = dir_figures,
        plot = gg_tmp + theme_paper + themes,
        width = width,
        height = height,
        device = cairo_pdf)
      # extrafont::embed_fonts(paste0(dir_figures, name, ".pdf"))
  }

# Load electricity rate data ---------------------------------------------------
  # Load the utility-specific datasets
  rt_pge <- dir_csv %>% paste0("electricityRatesPGE.csv") %>%
    read_csv() %>% data.table()
  rt_socal <- dir_csv %>% paste0("electricityRatesSocal.csv") %>%
    read_csv() %>% data.table()
  # Grab desired variables from SoCal's data
  rt_socal <- rt_socal[, c(1:2, 19:23, 18)]
  # Convert dates
  rt_pge[, `:=`(
    date_start = mdy(date_start),
    date_stop = mdy(date_stop)
  )]
  rt_socal[, `:=`(
    date_start = mdy(date_start),
    date_stop = mdy(date_stop)
  )]
  # Bind together
  rt_dt <- rbindlist(list(rt_pge, rt_socal), fill = T, use.names = T)
  # Count the number of days in each rate schedule
  rt_dt[, days := (date_stop - date_start + 1) %>% as.numeric()]
  # Expand to day level
  rt_dt <- rt_dt[rep(seq.int(1, nrow(rt_dt)), rt_dt$days), 1:9]
  # Add day-of-rate
  rt_dt[, rate_day := seq(1, .N), by = .(utility, date_start)]
  # Add day variable
  rt_dt[, date := date_start + rate_day - 1]
  # Add month, year, and month-year
  rt_dt[, `:=`(
    mo = month(date),
    yr = year(date)
  )][, ym := paste0(yr, mo)]
  # Keep first of the month
  el_dt <- rt_dt[day(date) == 1, .(month_date = date, utility, tier1, tier3)]
  # Stack
  el_dt <- rbindlist(list(
    el_dt[, .(month_date, utility, price = tier1, price_type = "elec_1")],
    el_dt[, .(month_date, utility, price = tier3, price_type = "elec_3")]
  ), use.names = T, fill = T)

# Load gas rate data -----------------------------------------------------------
  # Load data
  gas_dt <- readRDS(paste0(dir_fig, "ggplotPriceSeries.rds"))$data
  # Wide to long
  gas_dt <- gas_dt[, .(month_date, charge_base, charge_excess, utility)] %>%
    melt(id.var = c("month_date", "utility"),
    variable.name = "price_type", value.name = "price")
  # Keep the only the base rate
  gas_dt <- gas_dt[price_type == "charge_base"]
  # Create joint variable of utility and price type
  gas_dt[, price_type := "gas"]

# Join gas and electricity data ------------------------------------------------
  # Merge datasets for regression
  reg_dt <- merge(
    x = el_dt, all.x = T,
    y = gas_dt[, -"price_type"], all.y = F,
    by = c("month_date", "utility"),
    suffixes = c("_e", "_g"))
  # Stack datasets for plotting
  price_dt <- rbindlist(list(gas_dt, el_dt), use.names = T, fill = T)
  # Limit to 2010 through 2014
  price_dt <- price_dt[year(month_date) %>% between(2010, 2014)]
  # Add month, year, and month-year
  reg_dt[, `:=`(
    mo = month(month_date),
    yr = year(month_date)
  )][, ym := paste0(yr, mo)]
  price_dt[, `:=`(
    mo = month(month_date),
    yr = year(month_date)
  )][, ym := paste0(yr, mo)]
  # Create utility by price type
  price_dt[, utility_price_type := paste0(utility, "_", price_type)]
  # Create utlity comparison dataset for plotting differences
  wide_dt <- dcast(price_dt,
    month_date ~ utility + price_type, value.var = "price")
  # Take within-month, cross-utility differences (pge - socal)
  wide_dt[, `:=`(
    gas_diff = pge_gas - socal_gas,
    elec1_diff = pge_elec_1 - socal_elec_1,
    elec3_diff = pge_elec_3 - socal_elec_3
  )]
  # To long for a differenced data table
  diff_dt <- wide_dt[, c(1, 8:10)] %>%
    melt(id.var = "month_date", variable.name = "type", value.name = "diff")
  # Demean
  diff_dt[, diff_dm := diff - mean(diff), by = type]

# Plot: Times series of levels -------------------------------------------------
  gg_tmp <- ggplot(data = price_dt,
    aes(x = month_date, y = price,
      linetype = utility_price_type,
      color = utility_price_type)) +
    geom_line(size = 0.35) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = min(price_dt$month_date),
      color = "grey20", size = 0.4) +
    # theme_beamer +
    # ggtitle("Time series of electricity and gas prices for PG&E and SoCal",
    #   subtitle = "Tiers 1 and 3 for electricity; tier 1 for natural gas") +
    xlab("Date (month)") +
    ylab("Price") +
    scale_y_continuous(
      limits = c(0, 1.3),
      labels = scales::dollar
    ) +
    # scale_color_manual("Utility and type of price:",
    scale_color_manual("Price:\nUtility, type, tier",
      values = rep(c("#6A1B9A", "#ffc107"), each = 3),
      labels = c("PG&E: Elect. Tier 1", "PG&E: Elect. Tier 3", "PG&E: Gas Tier 1",
        "SoCal: Elect. Tier 1", "SoCal: Elect. Tier 3", "SoCal: Gas Tier 1")) +
    # scale_linetype_manual("Utility and type of price:",
    scale_linetype_manual("Price:\nUtility, type, tier",
      values = rep(c("solid", "dotted", "twodash"), 2),
      labels = c("PG&E: Elect. Tier 1", "PG&E: Elect. Tier 3", "PG&E: Gas Tier 1",
        "SoCal: Elect. Tier 1", "SoCal: Elect. Tier 3", "SoCal: Gas Tier 1")) +
    guides(col = guide_legend(byrow = T, nrow = 2))

  # save_beamer(gg_tmp, "gasElecSeriesLevels",
  #   width = 16 * 0.85, height = 10 * 0.85,
  #   themes = theme(
  #     legend.direction = "horizontal",
  #     legend.title = element_text(face = "bold"),
  #     legend.text = element_text(size = 15),
  #     legend.key.width = unit(1.5, "cm"),
  #     legend.key.height = unit(0.75, "cm"),
  #     NULL))

  # # Save for paper: Full size
  # save_paper(gg_tmp, name = "gasElecSeriesLevels",
  #   width = 9, height = 5.6,
  #   themes = theme(
  #     legend.direction = "horizontal",
  #     legend.key.height = unit(0.75, "cm"),
  #     NULL)
  # )

  # Save for paper: Small size
  save_paper(gg_tmp, name = "gasElecSeriesLevelsSmall",
    width = 6.4, height = 3.8,
    themes = theme(
      legend.direction = "horizontal",
      legend.key.width = unit(0.75, "cm"),
      legend.key.height = unit(0.5, "cm"),
      legend.title = element_text(size = 10),
      legend.text = element_text(size = 10),
      # legend.margin = margin(t=-0.75, r=0.1, b=-0.5, l=0.1, unit="cm"),
      NULL)
  )

# Plot: Times series of differences --------------------------------------------
  # Time series of differences
  gg_tmp <- ggplot(
    # data = diff_dt,
    data = diff_dt[type != "elec3_diff"],
    aes(x = month_date, y = diff_dm,
      linetype = type,
      color = type)) +
    geom_line(size = 0.35) +
    geom_hline(yintercept = 0, color = "grey20", size = 0.4) +
    geom_vline(xintercept = min(price_dt$month_date),
      color = "grey20", size = 0.4) +
    # theme_beamer +
    # ggtitle("Time series of the within-month differences between PG&E and SoCal",
    #   subtitle = "Difference in first-tier prices for electricity and gas (PG&E minus SoCal)") +
    xlab("Date (month)") +
    ylab("Price difference") +
    scale_y_continuous(
      labels = scales::dollar
    ) +
    scale_color_manual("Industry difference (PG&E price minus SoCal price):",
      # values = c("#FC845F", "#6A1B9A"),
      values = c("grey75", "grey40"),
      labels = c("Gas", "Electricity")) +
    scale_linetype_manual("Industry difference (PG&E price minus SoCal price):",
      values = c("solid", "twodash"),
      labels = c("Gas", "Electricity"))

  # save_beamer(gg_tmp, "gasElecSeriesDiff",
  #   width = 16 * 0.85, height = 10 * 0.85,
  #   themes = theme(
  #     legend.direction = "horizontal",
  #     legend.title = element_text(face = "bold"),
  #     legend.text = element_text(size = 15),
  #     legend.key.width = unit(1.5, "cm"),
  #     legend.key.height = unit(0.75, "cm"),
  #     NULL))

  # Save for paper: Full size
  save_paper(gg_tmp, name = "gasElecSeriesDiff",
    width = 9, height = 5.6,
    themes = theme(
      legend.direction = "horizontal",
      legend.key.height = unit(0.75, "cm"),
      NULL)
  )

  # Save for paper: Small size
  save_paper(gg_tmp, name = "gasElecSeriesDiffSmall",
    width = 6.4, height = 3.8,
    themes = theme(
      legend.direction = "horizontal",
      legend.key.width = unit(0.75, "cm"),
      legend.key.height = unit(0.5, "cm"),
      legend.title = element_text(size = 10),
      legend.text = element_text(size = 10),
      # legend.margin = margin(t=-0.75, r=0.1, b=-0.5, l=0.1, unit="cm"),
      NULL)
  )

# Regression -------------------------------------------------------------------
  #
  felm(price_e ~ price_g,
    data = reg_dt[price_type != "elec_3" & utility == "pge"]) %>% summary()
  felm(price_e ~ price_g,
    data = reg_dt[price_type != "elec_3" & utility == "socal"]) %>% summary()
  felm(price_e ~ price_g | utility + ym,
    data = reg_dt[price_type != "elec_3"]) %>% summary()
  #
  felm(price_e ~ price_g,
    data = reg_dt[price_type != "elec_1" & utility == "pge"]) %>% summary()
  felm(price_e ~ price_g,
    data = reg_dt[price_type != "elec_1" & utility == "socal"]) %>% summary()
  felm(price_e ~ price_g | utility + ym,
    data = reg_dt[price_type != "elec_1"]) %>% summary()
